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Abstract 

The structural and elastic properties of orthorhombic black phosphorus have been investigated 
using first-principles calculations based on density functional theory. The structural parameters 
have been calculated using the local density approximation (LDA), the generalized gradient ap- 
proximation (GGA), and with several dispersion corrections to include van der Waals interactions. 
It is found that the dispersion corrections improve the lattice parameters over LDA and GGA in 
comparison with experimental results. The calculations reproduce well the experimental trends un- 
der pressure and show that van der Waals interactions are most important for the crystallographic 
6-axis, in the sense that they have the largest effect on the bonding between the phosphorus layers. 
The elastic constants are calculated and are found to be in good agreement with experimental 
values. The calculated C22 elastic constant is significantly larger than the Cn and C33 parameters, 
implying that black phosphorus is stiffer against strain along the a-axis than along the b- and 
c-axes. From the calculated elastic constants, the mechanical properties such as bulk modulus, 
shear modulus, Young's modulus and Poisson's ratio are obtained. The calculated Raman active 
optical phonon frequencies and their pressure variations are in excellent agreement with available 
experimental results. 
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I. INTRODUCTION 



Phosphorus, an extraordinary element in the group 5A of the Periodic Table, exists 
in different allotropic modifications including white, black and amorphous red at ambient 
conditions.- Among these allotropic forms black phosphorus is the most stable one.- It is 
a semiconductor- that crystallizes in an orthorhombic structure.- At ambient conditions, 
orthorhombic black phosphorus has a puckered layered structure where each phosphorus 
atom is covalently bonded with three neighboring phosphorus atoms in the same layer, the 
layers being weakly coupled by van der Waals (vdW) forces. At high pressures orthorhombic 
black phosphorus transforms to a rhombohedral structure around 5 GPa, and transforms 
further on to a simple cubic structure at 10 GPa.-~- At very high pressures of 137 GPa and 
262 GPa, further structural phase transitions have been studied by Akahama et al.— ^ using 
x-ray diffraction. The measurements of the physical properties of black phosphorus such as 
the anisotropic electrical resistivity- and the optical reflectance have been reported,— 1 ^ as 
have been angle resolved ultraviolet photoemission spectra,— infrared absorption spectra,— 
and phonon dispersions obtained by inelastic neutron scattering. 16 

On the theoretical side, investigations of the electronic structure^"— using the linear 
muffin tin orbital (LMTO) method,— and optical properties using the self consistent pseu- 
dopotential method were published.-^ Lattice dynamical properties of black phosphorus 
including phonon dispersion relations were obtained™ using a valence force field method. 
A number of studies discuss the high-pressure behavior of black phosphorus, notably the 
structural stability and the electronic structure.—^— Recently, Du et al.— performed cal- 
culations for single layered and bulk black phosphorus and reported that black phosphorus 
may be used as anode material for lithium ion batteries. 

Since black phosphorus is a layered structure material that binds through van der Waals 
interactions between the layers it is necessary to include the vdW interactions in the theo- 
retical calculations of the ground state structural properties of this system. Van der Waals 
forces, or dispersive forces, result from the interaction between fluctuating multipoles with- 
out the need of an overlap of the electron densities. The standard density functional theory 
(DFT) functionals such as the local density approximation (LDA), the generalized gradient 
approximation (GGA), and various hybrid functionals do not describe properly the vdW 
forces, which are non-local in nature. To overcome this difficulty, we have used dispersion 
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corrected density-functional theory (DFT-D) approaches, which enhance the total-energy 
functional by a damped^ 3 - interatomic pairwise potential energy of the form CeR" 6 , where 
R denotes the interatomic distance. The main focus of the present work is to calculate 
the ground state structural and dynamical properties of black phosphorus, and to compare 
several schemes for inclusion of vdW interactions. 

The paper is organized as follows: The next section presents the computational methods 
used, and the results and discussions are presented in section III. The last section, IV, 
contains a summary. 



II. COMPUTATIONAL DETAILS 

The calculations were performed using the CAmbridge Series of Total Energy Package 
(CASTEP) 34 ^ code based on DFT using Vanderbilt-type ultra soft pseudo potentials 36 for 
the electron-ion interactions. The exchange-correlational potential of Ceperley and Alder— 
as parameterized by Perdew and Zunger— (CAPZ) in LDA, and the GGA schemes of Perdew 
and Wang22 (PW91) and of Perdew-Burke-Ernzerhof^, (PBE) were used. Different DFT- 
D approaches to treat vdW interactions were employed, notably the Ortmann, Bechstedt 
and Schmidt 41 (OBS) correction to PW91, as well as the Tkatchenko and Schemer 42 (TS) 
and Grimme 3 ^ corrections to PBE. In the three methods, a correction of the form = 
f(Rij) x x R~ & is added to the DFT total energy for each pair (ij) of atoms separated by 
a distance Rij. f(R) is a damping function which is necessary to avoid divergence for small 
values of R, and Cq is the dispersion coefficient for the atom pair (ij). The three methods 
use a slightly different damping function (see Refs. |4l|, |33j and |42j for details), but while 
the Cq coefficients in the Grimme and OBS methods depend only on the chemical species, 
the TS method calculate them for every structure by using the fact that the polarisability 



depends on the volume (see Ref. |42|). In the present study, we did not attempt to tweak the 
various coefficients of the different methods and kept them as in the original publications. 
In particular, in black phosphorus enters only one Cq coefficient, for the P-P interaction, 
with values of Cq = 103.6, 81.26, and 110.5 eV-A 6 in the OBS, Grimme and TS (unsealed) 
approaches, respectively. These methods have shown to be realiable and are often used for 
electronic structure calculations. For instance, Grimme's method has been used recently to 
study the binding properties of various solids 4 ^ and the structure and dynamics of molecular 



spin crossover compounds.— The OBS method is also widely used, as exemplified by the 
study— of the interaction between NO2 molecules and the Au-111 surface. The method of 
Tkatchenko and Schefner is also very popular and was used recently to study the hydrogen 
bonds of ice under pressure^ and the interlayer sliding energy landscape of hexagonal boron 
nitride.— 



To perform our calculations, a plane wave cut-off energy of 520 eV is used and the first 
Brillouin zone of the unit cell is sampled according to the Monkhorst-Pack scheme,— by 
means of a (9x6x8) k-point sampling. The Pulay 4 ^ density mixing scheme was applied 
for the electron energy minimization process. During the structural optimization process 
iterations were continued until the change in total energy was less than 5 • 10~ 6 eV/atom 
and the maximum displacement of atoms less than 5 ■ 10 -5 A. The elastic constants are 
calculated using the volume-conserving strain technique as implemented in the CASTEP 
code.— 



At ambient conditions black phosphorus crystallizes in a base centered orthorhombic 
structure (see Fig. 1) with the space group of Cmca (no. 64). The experimental lattice 
parameters 4 - are a = 3.31 A, b = 10.47 A and c = 4.37 A. The crystal structure consists 
of puckered layers of P atoms with three short covalent bonds (two at 2.22 A, and one at 
2.28 A), with vdW bonding along the 6-axis, the vertical distance between layers being 3.27 
A. The orthorhombic unit cell contains 8 phosphorus atoms in the crystallographic (8/) 
positions 4 (0,u,v). The nearest and next-nearest atomic separations between the layers (c?i 
and g?2 in Fig. 1) are 3.6 A and 3.8 A. 



In the orthorhombic crystal structure there are nine independent elastic constants^: Cn, 
C22! C33, C44, C 55 , C 66 , C12, C13, C23. Here, we follow the conventions of the experimental 
works^£ in which the indices 1, 2, 3 are defined to correspond to the crystallographic direc- 
tions X, Y, Z (c, a, b), respectively, as indicated in Fig. 1. To calculate the elastic constants, 
the equilibrium crystal structure is deformed by the appropriate strain and the variation in 
total en ergy gives the elastic constants. This is explained in detail for orthorhombic crystals 



in Ref. 



51 



The phosphorus internal coordinates u and v are freely varied during these 



deformations. 



III. RESULTS AND DISCUSSIONS 



A. Structural properties 



The theoretical equilibrium crystal structure was obtained by full structural optimization 
including lattice constants and the internal coordinates u and v, using LDA, GGA, and the 
three considered versions of GGA-D. The calculated structural parameters are presented in 
Table I. Taking the experimental ambient volume as reference, we find the difference between 
calculated and experimental volume of -11 % with LDA, +8 % with PW91, and +10 % with 
PBE. The large errors reflect that none of these approximations to the exchange-correlation 
potential are able to capture correctly the nature of the interactions in black phosphorus. 
The DFT-D schemes perform much better (see Table I), the errors in predicted equilibrium 
volumes being only +1.5 %, +2.9 %, and -0.3 %, respectively, for the PW91-OBS, PBE-TS, 
and PBE-Grimme parametrizations. Considering the three lattice parameters, it is clear that 
the largest error is found for the 6-parameter, the characteristic dimension in the Z-direction, 
i.e. perpendicular to the covalent P layers. Thus, the interlayer interaction is dominated by 
dispersive forces (vdW), and an accurate description of this type of interaction is important. 
Comparison of the data for GGA-PBE and GGA-PBE-Grimme in Table I shows that the 
inclusion of vdW interaction reduces the c-parameter by ~ 3 %, while the 6-parameter 
is reduced by ~ 7 %. On the other hand, the internal parameters u and v are about 
equally well predicted by all functionals considered, probably because these parameters are 
determined by the P-P covalent bonds, which are usually well described by LDA and GGA. 
The overall best agreement with experimental parameters is achieved with the PBE-Grimme 
functional. Our LDA calculations are in reasonable agreement with the LDA calculations 
performed by Ahujap^ (see Table I). The minor differences are most likely due to the different 
parametrizations of the LDA functional used and to the different basis sets (plane waves in 
the present work versus full potential LMTO in Ref. 

To investigate the response of black phosphorus to external pressure, hydrostatic com- 
pression was applied to the unit cell in the pressure range of to 5 GPa in steps of 0.5 GPa. 
This was done by variable cell optimization under the constraint of a diagonal stress tensor 
with values of the diagonal elements fixed to specify the desired pressure. The change of 
the lattice parameters a, b and c with pressure is compared with the experimental values 55 



191). 




FIG. 1. (Color online) Crystal structure of black phosphorus. The nearest and next-nearest 
inter layer atomic distances are shown (d± and c^). 

and shown in Fig. [2j It is seen that the a lattice parameter shows only a weak variation 
under pressure, which reflects the strong bonding along this direction. This fact is correctly 
captured by all the functionals that we have used, although we obtain a larger variation 
than in the experiments. In fact all functionals predict a slight increase in a with pressure. 
The variation is much larger for b and c. For the b parameter, a reduction of ~ 0.52 A is 
observed experimentally over the pressure range of 0-4.5 GPa.— This is correctly reproduced 
when the vdW correction is included, but it is about twice as large with the pure PW91 and 
PBE functionals. As for the c parameter, experimentally a reduction by approximately 0.27 
A is seen in this pressure range (Fig. 2(c)). This variation is overestimated by all methods 
that we have tested, even those corrected for vdW interactions. 

These results are also quantified in terms of the pressure coefficients defined as (with x 
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TABLE I. The calculated ground state properties of orthorhombic black phosphorus at ambient 
pressure, a, b and c are the lattice parameters, V the volume of the orthorhombic unit cell, and u 



and v t 
37 and 



re internal crystallographic parameters. CAPZ denotes t 



re LDA parametrization of Refs. 



3gJ, PW91 and PBE the GGA parametrization of Refs 



39 and 



van der Waals parametrization, the parametrizations of Ref. 



41 



40, respectively. For the 



(OBS), Ref. 



42 



33 



(TS) and Ref. 

(Grimme) are considered. Also included are the results of the LDA calculations by Ahuja^ with 
the Hedin-Lundqvist parametrization (LDA-HL) 54 . 



Method 


XC 


a(A) 


b(A) 


c(A) 


v(A 3 ) 


u 


V 


DFT 


LDA-CAPZ 


3.28 


10.10 


4.07 


134.9 


0.1059 


0.0711 




LDA-HL a 


3.24 


10.19 


4.24 


140.0 


0.1013 


0.0789 




GGA-PW91 


3.29 


11.07 


4.52 


164.5 


0.0950 


0.0871 




GGA-PBE 


3.28 


11.22 


4.54 


167.1 


0.0935 


0.0876 


DFT-D 


GGA-PW91-OBS 


3.29 


10.63 


4.41 


154.1 


0.0996 


0.0832 




GGA-PBE-TS 


3.29 


10.82 


4.39 


156.3 


0.0979 


0.0830 




GGA-PBE-Grimme 


3.30 


10.43 


4.40 


151.3 


0.1017 


0.0833 


Exp ft 




3.3133 


10.473 


4.374 


151.77 


0.1034 


0.0806 


a : Ref. ] 





being either a, b, or c) 



b : Ref. 



7(X) = % 



(1) 



p=o 



which are shown in Table [TTl The pressure coefficient for the a lattice parameter is very 
small for all approximations, in accordance with the experimental data.- For the b lattice 
parameter, the pressure coefficient varies significantly between the different approximations, 
but comes very close to the experimental value with Grimme's corrected functional. The 
magnitude of the pressure coefficient for the c parameter is roughly twice as large as found 
experimentally, with the best agreement found again with Grimme's method. 

The calculated pressure- volume (PV) relation of black phosphorus is shown in Fig. |3]and 
compared to experimental results.— Apart from the offset in equilibrium volumes discussed 
above, it is noticable that both the GGA and the GGA-D calculations find black phosphorus 
more compressible than observed experimentally. 



The slopes of the calculated V(P) relations at P = are larger than seen in experiments,— 
with the exception of the Grimme approach. This implies that the bulk moduli calculated 
with these approaches are significantly smaller than the experimental value. For the Grimme 
functional the initial slope is very similar to the experimental value, and the calculated bulk 
modulus within the Grimme approach is B = 31.3 GPa, which is in good agreement with 
experimental values of 32.5 GPa 4 and 36 GPa.— Therefore, comparing all the theoreti- 
cal approaches applied here for black phosphorus, Grimme's functional shows overall best 
performance, and it is noticeable that the dispersion correction provided by the Grimme 
functional has the right magnitude at low pressure, but seems to be slightly too strong for 
elevated pressures. However, one should bear in mind that the experiments under pressure 
are performed at room temperature, while the theory pertains to T = 0. 

The variation under pressure of the internal coordinates u and v are illustrated in Fig. 
HI All functionals capture the increasing (decreasing) trends of u(P) (v(P)) found in the 
experiment.— Clearly, the LDA shows the poorest agreement with the experimental data, 
and the vdW corrected GGA show marginally better agreement with experiments than the 
pure GGA functionals. Also, Fig. [5] illustrates the pressure dependance of the nearest and 
next nearest interlayer atomic distances, d\ and di in Fig. [TJ It is seen that LDA grossly 
underestimates and the pure GGA functionals overestimate these parameters. For the vdW 
corrected functionals, the agreement with experiments is much better, and the pressure 
evolution of the inter-layer atomic distances is best reproduced with the Grimme functional. 
This is of course in line with the better agreement found for the equilibrium volume and the 
other properties under pressure. 



B. Elastic properties 

Elastic constants are fundamental mechanical parameters for crystalline materials and 
describe the stiffness of the material against externally applied strains. In the present 
study, we have calculated elastic properties for single crystal as well as polycrystalline black 
phosphorus within the GGA and vdW corrected GGA approximations. The calculated 
elastic constants are presented in Table III. All the calculated elastic constants satisfy Born's 



mechanical stability criteria 56 implying that the system is in a mechanical 
The elastic constants have been determined for black phosphorus in Refs. 



stable state, 
and 



52 



53J by 



TABLE II. The calculated first order pressure coefficients of the lattice parameters (Eq. (JTJ) ) of 
orthorhombic black phosphorus, using LDA-CAPZ as well as GGA with and without dispersive 
corrections (see caption of Table I). The experimental results are taken from Ref. |4j. Units are 
10 _3 A-GPa -1 . 



DFT 



DFT-D 



Expt. (Ref. 4) 



axis LDA PW91 PBE 



PW91-OBS PBE-TS PBE-Grimme 



b 
c 



11 

-123 
-101 



~ 1 
-460 
-125 



-566 
-135 



3 

-261 
-118 



2 

-411 
-147 



~ 
-124 
-82 



~ 
-125 
-55 



ultrasound velocity measurements. Following the convention of these works, the Cn, C22, 
and C33 elastic constants directly relate to sound propagation along the crystallographic c, a, 
and 6-axes, respectively, and reflect the stiffness to uniaxial strains along these directions.— 
The reported experimental values show some scatter, which probably reflects the difficulty 
of the measurements. Also, the elastic constants show some variation with temperature, as 



discussed in Ref. 



53 



From Table III, it is seen that the calculated elastic constants show some variations when 
evaluated with different functionals, which is a combined effect of the different effective 
forces and the different equilibrium structures obtained with the various functionals. The 
C22 elastic constant is significantly larger than Cn and C33. This shows that the material is 
stiffer for strains along the a-axis than along the b- and c-axes, consistent with the results 
shown in Fig. [2J which again can be attributed to the strong covalent bonds between 
P atoms along the crystallographic a-axis. For C22, the calculated values are consistent 



between the different functionals and in excel 
Ref. I52I, while the experimental value of Ref. 



ent agreement with the experimental value of 
53| is significantly larger. The Cn is relatively 



constant among the considered functionals and again in fair agreement with the experiment 
of Ref. 



52, while the value obtained in Ref. 



53 



is larger. The C 33 describes stiffness to stress 



along the 6-axis, and this has a small value when evaluated in pure GGA, while the vdW- 
corrected methods give larger values, for the OBS and TS flavors in reasonable agreement 
with the experimental values, while the Grimme approach significantly overestimates this 
parameter. The C44 and C55 constants are relatively small and in good agreement between 

in 
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FIG. 2. (Color online) (a), (b) and (c): pressure dependence of lattice parameters a, b and c up 



to 5 GPa within LDA, GGA and GGA-D, compared to experimental values (Ref. 
different scales on the vertical axes, in particular the small scale in (a). 



55|). Note the 



theory and experiment (although the value of C55 calculated with PBE-TS stands out as 
too large). On the other hand, for Cqq the calculated values are fairly consistent and in 
good agreement with the experimental value of Ref. 



531 . while the measured value reported 



1 1 




P (GPa) 



FIG. 3. (Color online) Comparison of the pressure dependence of the crystal volume as calculated 
within LDA, GGA and GGA-D, compared to experimental volumes (Ref. 



by Ref. |52| is a factor 3-4 smaller. Among the off-diagonal elastic constants, C12 is large 
and relatively consistent between the various functionals, while the C13 and C23 are small, 
and in some cases take negative values. There are no reported experimental values of these 
quantities. From the single crystal elastic constants the single- crystal bulk modulus B may 
be found.— This gives a value of B = 30.7 GPa (using the Grimme functional), which is in 
good agreement with the experimental values of B = 32.5 GPa£ or B = 36 GPa,— as well as 
with the value derived from the V(P) relation in Fig. [3l B = 31.3 GPa (also in the Grimme 
approximation). Altogether, the agreement between the measured and the calculated elastic 
constants is fair while not perfect. 

The bulk and shear moduli of p oh/crystalline black phosphorus may be obtained from 
the elastic constants using the Voigt-Reuss-Hill approximations.— ~— Here, the Voigt and 
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FIG. 4. (Color online) Internal coordinates, u (a) and v (b), of black phosphorus up to 5 GPa as 



calculated within LDA, GGA and GGA-D, compared to experimental values (Ref. 



3a). 



Reuss approximations represent extreme values, and Hill recommended that their average is 



used in practice 
given by Ref. 



51 



or polycrystalline samples. Their definitions for orthorhombic systems are 
and their values for black phosphorus are quoted in Table IV. Using the 



calculated bulk and shear moduli in the Hill approximation, Bu and Gh, we have evaluated 
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FIG. 5. (Color online) Calculated (a): nearest neighbour (di) and (b): next-nearest nieghbor (o^) 
interlayer atomic separations (cf. Fig. 1) of black phosphorus up to 5 GPa with and without 



dispersion correction and compared with experiments (Ref. 



55|). 



the Young's modulus, E, and Poisson's ratio, a. Young's modulus reflects the resistance of a 
material towards uniaxial tensions, while Poisson's ratio gives information about the stability 
against shear strain. The results obtained are summarized in Table IV. It is seen that for 
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polycrystalline black phosphorus, Bh > Gh, which implies that the parameter limiting the 
mechanical stability of black phosphorus is the shear modulus Gh- The bulk and shear 
moduli provides information about the brittle-ductile nature of materials through Pugh's 
ratio.— According to Pugh, a Bh/Gh ratio below (above) 1.75 indicates the brittle {ductile) 
nature of the material. For black phosphorus the calculated value is Bh/Gh = 1.31, i. 
e. indicating a brittle nature of the material. Generally, Poisson's ratio is connected with 
the volume change of a material during uniaxial deformation and the nature of interatomic 
forces. If a is 0.5, no volume change occurs, whereas lower than 0.5 means that large volume 
change is associated with elastic deformation. 51 In our calculation o = 0.30 is obtained, 
which shows that a considerable volume change can be associated with the deformation in 
black phosphorus. 

The anisotropy in bonding between atoms in different planes is related to an anisotropic 
thermal expansion. This anisotropy plays an important role for all orthorhombic systems, as 
they influence the durability of the material to repeated heating cycles. The shear anisotropy 
factors measure this effect in terms of the elastic constants (see definitions for orthorhombic 



crystals in Ref. |5l| and l6ll). If the anisotropy factor is unity it represents complete elastic 



isotropy, while values smaller or greater than unity measure the degree of elastic anisotropy. 
For black phosphorus the calculated shear anisotropy factors Ai is close to unity, which 
shows that the resistance to shears of (100) planes (i.e. planes having c as a normal) along 
the a and b directions are approximately equally. In contrast, the A 2 and A 3 anisotropy 
factors are very different from unity, implying that shears of planes orthogonal to a and of 
planes orthogonal to b both behave very anisotropically. 



C. Transition to the A7 structure 

Black phosphorus transforms under pressure from the orthorhombic structure (A17) to 
a rhombohedral structure (A7) at roughly P t = 5 GPa.-~- This transition was studied pre- 
viously by Schiferl^ using an empirical local pseudopotential, and was analysed by Burdett 
et al.— using the concept of orbital symmetry conservation. Chang and Cohen^ performed 
calculations using the pseudopotential method within the local density approximation. How- 
ever, they found it necessary to apply a shift (of 2.3 mRy per P atom) to their calculated 
total energy of the A7 phase in order to reproduce the experimental value of the transition 
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TABLE III. Single-crystal elastic constants (C^, in GPa) and bulk modulus (B in GPa) of black 
phosphorus. All quantities are calculated at the respective theoretical equilibrium volumes using 
the various GGA and GGA-D functionals. 





DFT 




DFT-D 






Expt. 


parameter 


PW91 


PBE 


PW91-OBS 


PBE-TS 


PBE-Grimme 


Exp 


Cn 


45.2 


40.7 


44.4 


36.8 


52.3 


55. l a , 80 b 


C22 


186.7 


183.4 


192.4 


185.9 


191.9 


178.6 a , 284 b 


C33 


19.2 


13.0 


40.5 


30.6 


73.0 


53.6 a , 57 b 


C44 


10.3 


8.6 


15.9 


31.2 


25.5 


21. 3 a , 11. l a , 17.2 b 


C55 


4.2 


2.6 


5.6 


23.4 


8.8 


5.5 a , 10.8 6 


C66 


48.3 


48.5 


57.3 


57.4 


63.6 


14.5 a , 15. 6 a , 59.4 b 


C12 


32.3 


30.0 


35.6 


31.5 


40.8 




C13 


-4.3 


-4.6 


-5.3 


-0.9 


0.6 




C23 


1.0 


-1.6 


1.5 


-0.6 


8.3 




B 


11.6 


8.1 


18.6 


16.2 


30.7 


32.5 C , 36 d 



a: Ref. |52| (room temperature); b: As read from figures of Ref. l53l (T=0 K); c: Ref. |4|; d: From 



PF-data of Ref. 



55 



TABLE IV. Calculated values of the polycrystalline bulk (Bx) and shear moduli (Gx) in the Voigt, 
Reuss and Hill approximations {X = V,R,H, respectively). Units are GPa. Also listed are the 
Young's modulus, E (in GPa), Poisson's ratio, a, and the shear anisotropy factors, A». All values 
are calculated at the theoretical equilibrium volume using the PBE-Grimme functional. 



Parameter By 


Br 


B H 


Gy 




Gh 


E 


a Ai 


A 2 


A 3 


46.3 


30.7 


38.5 


37.4 


21.4 


29.4 


70.3 


0.30 0.8 


0.1 


1.6 



pressure, which is pointing out the inadequacy of the LDA to describe correctly the energet- 
ics of black phosphorus. This is in contradiction with the results of Ahujar^ who obtained 
a transition pressure of Pt= 4.5 GPa using the LMTO method and LDA. In the present 
work we also studied the transition pressure using Grimme's vdW functional. This leads to 
a value of P t = 0.8 GPa, implying that although the Grimme functional is correctly ranking 
the A17 structure versus the A7 structure in terms of the total energy at P = GPa, the 
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transition happens at a too low pressure in the calculation compared to the experimental 
pressure. Similar calculations were carried out using either of the PW91, PBE, OBS or TS 
functionals, but in all cases the transition pressure is found to be too low, within a 1 GPa 
range around the value obtained with the Grimme functional. 

Possibly, a better agreement with the experimental transition pressure could be obtained 
using a more advanced description of the van der Waals interactions, such as the adiabatic- 
connection fluctuation-dissipation theorem^— in the random phase approximation or quan- 
tum Monte-Carlo^ method. However, these methods are numerically heavier than the sim- 
ple dipole-dipole corrections considered here, and to our knowledge the implementation of 
the corresponding forces necessary to perform the structural relaxation under pressure is 
not available. 

D. Optical Phonons 

The optical phonons of black phosphorus have been obtained by the linear response 
method within density functional perturbation theory— In this method the force-constant 
matrix is obtained by differentiation of the Hellmann-Feynman forces on atoms with respect 
to the ionic coordinates. In order to understand the optical phonons of black phosphorus, it 
is necessary to first consider the dynamical properties of this structure. The primitive cell 
contains four atoms leading to 12 vibrational modes, i. e. 9 optical and 3 acoustic phonons 
branches. From the group theoretical analysis of the Cmca space group, the irreducible 
representation at the Brillouin zone center is 

^acoustic B\ u -\- B 2 u ~\~ B^ u 

and 

J- optical B\u + B 2u + A u + 2A g + Bi g + B 2g + 2B 3g 

All the optical phonons are Raman active except the B\ u and B 2u modes, which are 
infrared active, and the A u mode, which is silent. The optical phonons of black phosphorus 
were investigated by Sugai and Shirotani^ and by Akahama et al.— The Raman frequencies 
turn out to be relatively insensitive to the details of the functional employed, for which 
reason we show the results obtained with the GGA-Grimme fuctional. This is probably a 
reflection of the fact that the optical phonons are described by near-neighbour interatomic 

17 



forces, which are dominated by covalent interactions, which are generally well accounted 
for by the LDA. The eigenvectors of the Raman modes, as calculated at ambient pressure 
within the GGA-Grimme approach are illustrated in Fig. El in good agreement with earlier 
work.— ^ From the calculated Raman eigen-modes, the A 1 , A 2 g and B\ modes are seen 
to constitute inter-layer vibrations with a strong component along the b direction. In the 
modes B\ g and B2 9 the atoms move along the a direction, while in the B\ g mode the atoms 
move along the c direction. Ref. |26| studied the pressure variation of the Raman shifts. In 



Fig. [7] the calculated Raman frequencies under pressure up to 5 GPa are compared with 
these experimental results. The agreement between theory and experiment is excellent. The 
strongest pressure dependence is seen for the inter-layer A l g mode, for which the pressure 
coefficient is du/dP = 8.2 cm _1 /GPa in theory, and dui/dP = 5.78 cm _1 /GPa according to 
experiment.™ 



IV. CONCLUSIONS 

First-principles calculations based on the DFT have been used to investigate the struc- 
tural properties of black phosphorus within LDA and GGA. It was found that the calculated 
material parameters differ significantly from the experimental data, while including empirical 
methods to account for van der Waals interactions (DFT-D) leads to structural properties in 
good agreement with experiments. This demonstrates that van der Waals interactions play 
a major role in black phosphorus. Among the OBS, TS and Grimme parametrizations of van 
der Waals interactions, the calculated structural parameters are in closest agreement with 
experimental data when PBE-Grimme is used. Also, the influence of external pressure on the 
lattice parameters were calculated. The b lattice parameter was found to be most sensitive 
to which vdW functional was used, reflecting the dominance of van der Waals interactions 
along the crystallographic 6-direction. Further, the elastic properties of single crystal as 
well as polycrystalline black phosphorus were examined using the GGA-PBE-Grimme func- 
tional. Our calculations confirm that black phosphorus is a mechanically stable anisotropic 
material, and that the system is stiff est along the a-axis. The calculated shear anisotropy 
factors for different directions indicate strong anisotropies for shears perpendicular to the a 
and b axes. The polycrystalline bulk, shear and Young's moduli, as well as Poisson's ratio 
were calculated. From the ratio of bulk to shear modulus, black phosphorus is found to be 
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FIG. 6. (Color online) Snapshots of Raman active modes for primitive black phosphorus. The 
listed frequencies are as calculated within GGA-Grimme at P=0. 



brittle. Finally, the frequencies of Raman active vibrational modes were calculated and they 
are in excellent agreement with the available experimental data. 
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FIG. 7. (Colour online) Calculated (PBE-Grimme) Raman frequencies under pressure compared 
to experiments (Ref. bd ). 
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